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Abstract. The stability and convergence rate of Olver's collocation method for 
the numerical solution of Riemann-Hilbert problems (RHPs) is known to depend 
very sensitively on the particular choice of contours used as data of the RHR By 
C manually performing contour deformations that proved to be successful in the 

asymptotic analysis of RHPs, such as the method of nonlinear steepest descent, 
the numerical method can basically be preconditioned, making it asymptotically 

^3 stable. In this paper, however, we will show that most of these preconditioning 

deformations, including lensing, can be addressed in an automatic, completely 

algorithmic fashion that would turn the numerical method into a black-box solver. 

To this end, the preconditioning of RHPs is recast as a discrete, graph-based 
optimization problem: the deformed contours are obtained as a system of shortest 
paths within a planar graph weighted by the relative strength of the jump matrices. 
• The algorithm is illustrated for the RHP representing the Painleve II transcendents. 

t 

1. Introduction 

^ Remarkably many integrable problems in mathematics, mathematical physics, 

\Q and applied mathematics can be cast as Riemann-Hilbert problems (RHPs): classi- 

cal orthogonal polynomials and special functions, Painleve transcendents, nonlin- 
ear PDEs related to the inverse scattering transform, and distributions in random 
matrix theory as well as in random combinatorial problems, to name just a few |j5]. 
A fruitful point of view is that RHPs generalize the representation of classical 
^ special functions by contour integrals; like these they have extensively been used in 

i-H establishing deep asymptotic results such as connection formulae for the Painleve 

KjjJ transcendents [4]. Here, a fundamental tool is the method of nonlinear steepest de- 

• 1^ scent, which was introduced by Deift and Zhou [2] to the asymptotic analysis of 

/\ oscillatory RHPs. 

^ Only quite recently, starting with a novel direct spectral collocation method of 

Olver |6j, RHPs have become the subject of study in numerical analysis. It turned 
out that the stability of this numerical method and its approximation properties 
strongly depends on matching most of the steps from the asymptotic analysis of 
the problem at hand: "One can expect that whenever the method of nonlinear 
steepest descent produces an asymptotic formula, the numerical method can be 
made asymptotically stable'' [7, p. 2]. This way, the method is kind of hybrid: only 
after manually performing a series of expert analytic steps to deform the given 
RHP into another, equivalent, one, the thus "preconditioned'' RHP is taken as 
input to the numerical algorithm. Hence, the use of this numerical method has 
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been limited so far to an audience that would have competent operational access 
to such a kind of expert knowledge. 

In this paper we will give a proof of concept that most, if not all, of these defor- 
mations (at least if they were meant to stabilize the numerical method) can be 
addressed with an automatic, completely algorithmic approach that would turn 
the numerical method into a black-box solver for the userQ Following Bornemann 
and Wechslberger [i], who dealt with similar problems for contour integrals, we 
will recast the preconditioning of RHPs as a discrete, graph-based optimization 
problem: the desired deformation corresponds to a system of shortest paths within 
a weighted planar graphQ Though we are not able, at this stage of our study, to 
determine the precise complexity class of this particular discrete optimization 
problem (NP-hard, polynomial, etc.) or to prove that our polynomial greedy algo- 
rithm would approximate the optimum within a certain range (which would be 
sufficient for the purpose of preconditioning), we will demonstrate for the example 
of the Painleve II transcendents that, first, it will improve the stability of the input 
RHP significantly by several orders of magnitude and that, second, the resulting 
deformations of the RHP closely match what people have obtained by applying 
the method of nonlinear steepest descent. 

Riemann-Hilbert problems. To fix the notation, we consider RHPs for a given 
oriented contour T, which is a finite union of simple smooth curves Fy (; = 1, . . . , fc) 
in C. By removing the finitely many points of self-intersection from F we obtain 
F^ C F. Given a matrix- valued jump function G : F^ ^ GL(m,C), the RHP 
determines a holomorphic function O : C \ F ^ GL(m,C) satisfying 

^+{z) = ^-{z)G{z) (zGF°), 0(00) =F 

Here, 0^(z) denotes the non-tangential limit of O(z') as ^ z from the positive 
(negative) side of the contour. Existence and uniqueness of a solution O can be 
shown under some appropriate smoothness and decay assumptions on the jump 
function G [jjj. To simplify the discussion of contour deformations, we assume 
that there are entire functions Gy : C ^ GL(m,C) that continue the jump data G 
given on the part Fy of the contour F: 

^IrOnFy = GylrOnFy- 

We consider the pairs (Fy, Gy) (; = 1, . . . ,fc) as the data of the RHP. Most often one 
is not interested in the full solution 0(z) of the RHP but only on some derived 
quantities at 00, e.g., the residue 

res 0(z) = lim z(I - 0(z)). 

z=oo z^oo 

Example: Painleve II. Throughout this paper we will illustrate our ideas for the 
RHP representing the Painleve II equation 

Uxx = XU+ lu^. 



^To begin with, in this paper we study plain contour deformations and lensing; other important 
concepts of deformations, such as ^-functions, will be subject of subsequent refinements of our work. 
walk in a graph is a sequence of adjacent vertices, a path is a simple (non self-intersecting) walk. 
^The second condition is meant to imply that O has a holomorphic continuation at 00. 
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Figure i. The six rays Fy of the RHP representing Painleve IL 



The general solution u{x) = u{x; 81,82) of this second-order ODE in the complex 
domain will depend on two independent complex parameters Si and which 
are fixed in the following setup of the RHP: with the six rays (see Fig. [i| 

Fy = {se^M2y-i)/6 ; s > 0} (; = 1,. . .,6), 
parameters Sy (; = 1, . . . , 6) interrelated by 

Si - S2 + S3 + S1S2S3 = 0, S4 = -Si, S5 = -S2, S6 

and the jump matrices 



\ Sye-^(-)^ 

^ 1 0' 

^Sye+^(^) 1^ 



; even. 



; odd. 



with the phase function 



(z) = + lixz 



the solution O of the RHP yields 

uix;8i,82) = -2 res O1 2(z) = 2 lim zOi 2(z). 

z=oo ' z^oo 

Note that the independent variable x of the ODE enters the RHP as a parameter 
of the phase function 6: the RHP (with independent variable z) amounts thus for a 
pointwi8e evaluation of the Painleve transcendent z^(x;si,S2). 

The numerical method of Giver. Olver [6J constructed his spectral collocation 
method by recasting RHPs as a particular kind of singular integral equation. Upon 
writing 

0(z) = J + CrLI(z) 

with the Cauchy transform of a matrix-valued function U :T ^ GL(m,C), namely 



the RHP becomes the linear operator equatiori 

AU{z) = U{z) - Cf !J(z) • (G(z) -I) = G(z) - I. 



(1) 



4ln the singular case Si — S2 — ±z, there is a one-parameter family of solutions depending on S3. 
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Here, LI(z) denotes the non-tangential limit of CyU{z') as ^ z from the 
positive (negative) side of the contour; there is the operator identity Cy — = I- 
The residue at oo becomes simply the integral 

res 0(2) = -1. / 

z=co 2m Jt 

Without going into details, in this paper it suffices to note that the n-point numerical 
approximation of (jij yields a finite-dimensional linear system 

where the ;th component of the solution Un is a matrix that approximates U{zj) at 
the collocation point zy G F (; = 1, . . . , n). The stability of the method is essentially 
described by the condition number 

% = ^(^n) = ll^n^ll ■ ll^^ll 

of this linear system: altogether, one would typically suffer a loss of log^g Kn signif- 
icant digits. Under an additional assumption, which can be checked a posteriori 
within the numerical method itself, Olver and Trogdon \y, Assumpt. 6.1 and 
Lemma 6.1] proved a bound of the formj^ 

Kn = 0(k(A)) 

in terms of the condition number k{A) = • ||A|| of the continuous operator 

A, with constants that are midly growing in the number of collocation points n. 
Here, the operator norm of A is obtained by acting on L^(r). Extending Un to all 
of r by interpolation, Olver and Trogdon ^ Eq. (6.1)] also state an error estimate 
of the form 

ll^-^n|lL2(r) ^ ^^^(^)^^^^"^ll^llH^(r) 
with some j6 > 0. Since, for jump matrices G that are piecewise restrictions of 
entire functions, k can be chosen arbitrarily large, one gets spectral accuracy 

Preconditioning of RHPs. Thus, stability and accuracy of the numerical method 
depend on which blows up in many problems of interest. For instance, the 
undef ormed version of the RHP for Painleve II with the contours in Fig. [1] has 
Kn ~ 2.2 • 10^ for Si = 1 and S2 = 2 and x = —10 (see also Fig.jzjvor varying x). 

Now, it is important to understand that k{A) is the condition number of the 
RHP for the restricted data (F, G) but not for the jump data Gy (; = 1, . . . , fc) which 
are obtained from analytic continuation. If the continued data are explicitly given, 
and are not themselves part of the computational problem]^ it should be possible 
to deform the RHP to an equivalent one with data (f , G) and 

We call such a deformation a preconditioning one. In fact, Olver and Trogdon ^ 
argued that preconditioning is possible whenever the method of nonlinear steepest 
descent produces an asymptotic formula; Fig. [4] shows a typical sequence of such 
manually constructed preconditioning deformations for the Painleve II RHP. 

^They employ the estimate || A|| < \/2(l + ||G — J||^co(p) ||Cp ||) in tlieir statements. 
^Analytic continuation corresponds to solving a Cauchy problem for the elliptic Cauchy-Riemann 
differential equations; it is, therefore, an ill-posed problem. 
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Figure 2. Comparison of the condition number of the original contour 
and of a deformed contour optimized by the greedy algorithm of ^ 
for the Painleve II RHP with (si,S2) = (1/2). The condition number of 
the deformed contour is roughly constant for all values of x while the 
condition number of the original contour grows exponentially fast for 
decreasing values of x. Note that condition numbers larger than 10^^ 
(dashed line) obstruct the computation of even a single accurate digit in 
machine arithmetic, indicating severe numerical instability. 



Though it seems to be difficult to extract a single governing principle for all 
the ingenious deformations that are used in the asymptotic analysis of RHPs, we 
base our algorithmic approach on the following simple observation: if there are 
no jumps in the RHP, that is if G = J, we have A = I and therefore k{A) = 1. 
By continuity, G ^ I in some sufficiently strong norm would certainly imply 
k{A) 1, such that a reasonably small ||G — J|| will probably yield a moderately 
sized condition number k{A). We conjecture that such an estimate can be cast in 
the form 

KiA)^<pi\\G-I\\„s,p(^r)) 

for some Sobolev W^'^-norm and some monotone function (p that is independent 
of (T, G). A good preconditioning strategy would then be to make ||G — ^||iA^s,p(p) 
as small as possible, we call it the relative strength of the jump matrix G. 

In the lack of any better understanding of the precise dependence of k{A) 
on the RHP data (T, G) we suggest to use ||G — JH^^ij^p^ as a measure of relative 
strength: optimizing it led to significant reductions of the condition number in all 
of our experiments. However, the deformation algorithm itself will just use that 
the measure d{T;G) can be written as an integral over T, namely in the form 

d(r;G) = J^d{G{z))d\z\ 

for some function d : GL(m,C) [0, 00), which we call the local weight. 



6 



GEORG WECHSLBERGER AND FOLKMAR BORNEMANN 





c. (G^,f^) = LensingDeformation(G^,f ^) 
K Pi 250 



d. (G^,f3) = LensingDeformation(G2 f2^ 
K ^ 140 



Figure 3. Application of SimpleDeformation and LensingDeforma- 
TiON of ^to the Painleve II RHP with si = 1, S2 = 2 and x = -10. On 
the top left is the original contour of this RHP and the other contours are 
deformed versions of it. All contours have been calculated on a 17 x 17 
grid. The color encodes the magnitude of || G(z) — I||f , with green = 10~^^, 
yellow = 1, and red = 10^. The blue dots indicate the origin z = and the 
stationary points of the phase function 6 at z = x/2. The reduction 

of the condition number from (a) to (d) corresponds to an accuracy gain 
of about six digits. 



Preconditioning as a discrete optimization problem. Since our objective is pre- 
conditioning, the relative strength of the jump matrices does not really have to be 
minimized over all equivalent deformations (f , G) of a given RHP (F, G). For all 
practical purposes it suffices to consider just a very coarse, finite set of possible 
contours, namely paths within a planar graph. 

The basic idea is as follows: first, we restrict the problem to a bounded region of 
the complex plane and embed the part of the contour F belonging to that region as 
paths into a coarse, grid-like planar graph g = {V, E) (see Fig.[3]a for an example 
of the Painleve II RHP: because of a super exponential decay as z ^ 00 along each 
of the rays, G — I is already a computer zero outside the indicated rectangle). 

Second, for each the analytic continuation Gy of the jump data on Fy turns the 
graph g into an edge-weighted graph gj by using the (edge) weights 

dj{e)^ jd(Gj(z))d\z\ (eeE). 




Figure 4. Some manual constructions for the Painleve II RHP taken 
from Olver and Trogdon \y, p. 20]. Left: Deformation along the paths of 
steepest descent; right: deformation after lensing. The contours bifurcate 
at the stationary points of the phase function 6. 



Last, we replace Fy (within the bounded domain) by the shortest path (with the 
same endpoints as Fy) with respect to gy subject to the following constraint: the 
thus deformed RHP must be equivalent to the original one. 

It is this latter constraint which adds to the algorithmic difficulty of the prob- 
lem: the Fy cannot be optimized independent of each other. We will address this 
problem by a greedy strategy: the largest contribution to the weight constraints the 
admissible paths of the second largest one and so on; this will be accomplished by 
modifying the underlying graphs gy in the corresponding order. 

Fig. [3]b shows the result of such an algorithmic deformation for the Painleve 
II RHP {si = 1, S2 = 2, X = —10): the condition number is reduced by about six 
orders of magnitude. Further improvement is possible by performing a "lensing'' 
deformation, that is, by introducing multiple edges based on a factorization of G 
(see ^3.3] ). The results of two such steps are shown in Fig.[3]c and d (more steps 
would not pay off). Though the improvement of the condition number is more 
modest in these two steps, it is instructive to compare the algorithmic contour in 
Fig. [3]d with the manual construction of Olver and Trogdon [7] shown in Fig. [4] 

Fig. [2] compares, for varying values of x, the condition number of the original 
contour with that of the deformed contour optimized by the greedy algorithm of 
^ a uniform stabilization by preconditioning is clearly visible. 

Outline of the paper. In ^we discuss the two admissible deformations of RHPs 
that will be considered in this paper: simple deformations of contours and lensing 
deformations based on factorizations of the jump matrix G. We address the 
question of how to match the topological constraints of such deformations in the 
planar graphs attached to each part Fy of the contour. In ^we give an in-depth 
description (with pseudo code) of the greedy algorithm that aims at optimizing 
these deformations. Important steps are illustrated for the Painleve II RHP. Further 
details of the implementation are discussed in ^ 

2. Admissible Deformations of Riemann-Hilbert Problems 

We briefly recall two of the deformations that can be applied to RHPs. A more 
detailed description can be found in [4]. 
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undeformed 



deformed 



f 




r 



Figure 5. A simple deformation. 



2.1. Simple Deformations. Fig. [5] shows an example of such a deformation. In 
general, simple deformations allow to continuously move a contour part in the 
complex plane (thereby covering a region Q) as long as the following conditions 
are satisfied: 

(i) f does not cross other parts of F, 

(ii) O does not contain any other contour parts, 

(iii) G has a holomorphic continuation in Q. 

Then, the deformed RHP in Fig. [5] is solved by the function 



Conditions (i)-(iii) can be mapped to graph-constrained deformations as follows: 
Condition (i) can be handled by splitting a graph as shown in Fig. [6] if a path p 
corresponding to a part of a contour is given, like the path highlighted in blue, we 
duplicate the vertices of p and change all edges on the right side of p so that they 
are connected to the newly created vertices but not to the vertices of p itself. This 
way no path in the graph can cross p anymore. We will use gipirPir • • •] to denote 

a graph g which has been split in this fashion along the paths pi,p2f 

Condition (ii) is difficult to be built into the structure of a graph a priori, but 
it is easy to check for it a posteriori: the circle composed by F/ and f / should not 
enclose an endpoint of another arc. If violated, the algorithm simply stops (this 
never happened in our experiments; dealing with such a situation would require 
to break the deformations into smaller pieces). 

Condition (iii) can be handled by removing those regions from the graph where 
G does not have a holomorphic continuation. 

2.2. Multiple Deformations and Factorization: Lensing. Fig. [t] shows an exam- 
ple of such a deformation. To initialize, several copies of a contour part are created 
at one and the same location, where each copy corresponds to a factor of a given 
multiplicative decomposition of the jump matrix G. We call these copies the factors 
of this part of the contour F. These factors are then moved around in the complex 
plane subject to conditions (i)-(iii) and, additionally, the following condition: 
(iv) the mutual orientation of the factors must be preserved. 

For example, in Fig. [t] the order of the decomposition G = LDU requires that 
the factor LI is to the left of the factor D and that the factor D is to the left of the 
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g g[p] 

Figure 6. Illustration of split graphs. On the left side is the original graph 
g with a path p (highlighted in blue) along which it is about to be split. 
The split graph g[p] is on the right side, with all vertices and edges that 
have been changed or created highlighted in red. For a clear visualization 
the split in the graph has been enlarged by moving the vertex positions; 
in the actual graphs used by the algorithm the duplicated vertices stay 
at exactly the same position. All graphs are undirected, the arrow at the 
end of the blue path just indicates its orientation. 



undeformed 




G = LDU 



lensing 




Figure 7. A lensing deformation. 



factor L. To preserve this orientation in our deformation algorithm, we calculate the 
shortest path for just one of the factors. For the other factors we use a modification 



of the shortest enclosing circle algorithm of Provan ||8l, see 63.3 



3. The Greedy Algorithm 

3.1. Notation. 

• de : weight of the edge e 

• Pi + P2 : path Pi joined with path P2 

• V : reversed path P 

• P [u, v] : subpath from vertex u to vertex v within the path P 

• sp(g, u, v) : shortest path from vertex w to in the weighted graph g 
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• int(W) : homological interior of a closed walk W, that is, all vertices v of 
the graph that have winding number ind(W, i?) = ±1 w.r.t. W 

• p-/p-^ : path on the left/right side of the split along p in g[p] (see Fig.|6]) 



Algorithm i Optimized Simple Deformation 
i: procedure SimpleDeformation(G, F) 



2: n = I r I > n is the number of contour parts 

3: P = > fixed new paths 

4: p = > candidates for new paths 

5: f = > already processed contour parts 

6: Q = (1, . . . , n) > unprocessed contour parts 

7: g = > graphs corresponding to contour parts 

^- S"^ = ^ initial not split graphs 

9: for all z G Q do 

10: gi = graph with edge weights de = \\Gi — 

11: 17- = vertex in gi nearest to left endpoint of 

12: ■ = vertex in gi nearest to right endpoint of F/ 

13: end for 

15: while Q 7^ do 
16: for all z G Q do 

17: = sp(gf, i^-, i^-) > shortest path from v\ to i^- in 

18: end for 

19: zV = argmaxd(p/) 

20: P4 = P4 

Q = Q\(zV) 

22: for all z E f do > try to improve the path if there is an intersection 

23: if P4 n P/ 7^ then 

24: P = ImproveSharedSubpath(g'^, G, P, z, zV) 

25: end if 

26: end for 

27: P = PU(zV) 

28: for all z G Q do 

30: end for 

31: end while 

32: (G, f ) = MapToRHP(G, P) > create new contour parts 

33: return (G,f) 



34: end procedure 



3.2. Optimized Simple Deformations. The idea of Algorithm [1] goes as follows: 
first (lines 9-13), for each of the contour parts Fy and the corresponding jump matri- 
ces Gj (which are assumed to have a holomorphic continuation to the rectangular 
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Figure 9. Calculate the shortest paths for each Gy separately, highlighted 
in blue (line ' ' ^ " ^ 



17 



in Algorith m The one with the largest total weight 
is shown in magenta (line The color encodes the magnitude of 
|Gy(z) - I||f, with green = 10^^, yellow = 1, and red=10^^. 



region supporting the grid), a separate weighted graph gj with edge weights 

de = d{Gj{z)) d\z\ 

J e 

is created, see Fig.|8] Second (1 ines 16-18), each Fy is replaced by a shortest path 
that shares the same endpoints, see Fig. [9] The thus separately optimized paths, 
however, will in general not satisfy condition (i) of ^ that is, they will cross each 
other. Therefore, some of the paths have to be modified to match this condition, 
which increases the corresponding weight. By keeping, third (lines 19-20), the path 
P of dominant total weight fixed, we restrict such modifications to the other parts 
that contribute less to the condition number. By splitting, fourth (lines 28-29), all 
graphs along P and repeating (lines 16-18) the calculation of the shortest paths in 



the split graphs, we come up with paths that do not cross P, see Fig. 10 
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Figure lo. Recalculate the shortest paths (line[i7|in Algorithm [i) i n the 
graphs split along the optimal path of largest weight from Fig.fgj here 
shown in white. The new shortest paths are highlighted in blue, the one 
of maximal weight in magenta. 





Improvement of the subpath shared by the fixed white path 
(Algorithm [2|. After this improvement. 



Figure ii 

and the magenta path of Fig. 

both paths are fixed and the graphs are split along their union (the white 
Y-shape contour). From here it should be clear, how Algorithm [i] arrives 
(for a finer grid) at the deformed contours shown in Fig.|3]b. 



This procedure is then repeated, fifth (line 15), until all paths are fixed and, 
hence, non-crossing. (In each round of this loop another path gets fixed.) Finally, 
sixth (line [32]), the algorithm constructs the deformed contour data from the just 
calculated set of paths. For paths and subpaths that do not share an edge with 
another path we simply use the path and the corresponding jump matrix as new 
contour data. Subpaths which occur in more than one path will be mapped to 
new contour data by performing an "inverse lensing'': the new jump matrix is 
calculated as the (properly ordered) product of all the jump matrices sharing the 
that subpath. For example, if the paths P/ and Py have a common subpath s and P/ 
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is to the left of Py, the procedure MapToRHP creates a new contour part s with the 
jump matrix GyG/. 

In a situation as shown in Fig. [loj where the new optimal paths share a subpath 
with the already fixed ones, further improvement is possible (lines 23-25) by 
optimizing the shared subpath with respect to the weight obtained from combining 
the corresponding jump matrices (that is, the just mentioned "inverse lensing''). 
This procedure (Algorithm |2]| is schematically illustrated in Fig. 
of this procedure to the example of Fig. 



10 



is shown in Fig. 
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the application 



Algorithm 2 Improve Paths with a Shared Subpath 



procedure lMPROVESHAREDSUBPATH(g, G, P, z'l, Z2) 

= gi^ > create a temporary graph 

if Pi^ left of P,2 then 

set weights for edges e of g' to d{e) = ||G/2G/^ ~ 
end if 

if P/^ right of P/2 then 

set weights for edges e of g^ to d{e) = HG/^G/^ — 
end if 

s = Pi^ n P/2 > split Pi^ and P/^ into left,common and right subpath 

g'^g'[p',pi+K'Ph+K^ 

Pi, = sp(g',si,s_i) > = shortest path between start and end vertex of s 
Pi^=PlUp,UPJ^ 

if containsCircle(P/J then 
Pi^ = dropCircle(P/J 
s = P- n P- 

P^,=%{g',[P\Pi,i{P^,)l.{Pi,)-^ 

if Pi^ n P/2 7^ s then 

P = ImproveSharedSubpath(g, G, P, z'l, Z2) 
end if 

else if containsCircle(P/2) then 
P/2 = dropCircle(P/2) 
s = P- n P- 

if Pi^ n P/2 7^ s then 

P = ImproveSharedSubpath(g, G, P, z'l, Z2) 
end if 
end if 
return P 
end procedure 



14 



GEORG WECHSLBERGER AND FOLKMAR BORNEMANN 



4 


8 


12 


i6 


20 


4 


8 


12 


16 


20 


3 


7 


11 


15 


19 


3 


7 


f < 

11 


15 


19 


2 


6 


10 


14 


18 


2 


6 


10 


14 


18 


1 


5 


9 


13 


17 


1 


5 


9 


13 


17 


4 


8 


12 


16 


20 


4 


8 


12 


16 


20 


3 


7 


11 


15 


19 


3 


7 


II 


15 


19 


2 




r°, 


14 


18 


2 


1 1 

6 


10 


14 


18 


1 


5 


9 


13 


17 


1 


5 


9 


13 


17 



Figure 12. An illustration of applying one round of ImproveShared- 
SuBPATH (Algorithm [2|. Top left: a graph in which the two paths Pj and 
Pj shown in blue and green share the subpath (7,11,15). Top right: after 
calculating a shortest path from 7 to 15 for the combined weight (lines 
3-8 in Algorithm [2J, the blue path contains the circle (6,7,6). Bottom left: 
this circle gets removed from the blue path. Bottom right: an updated 
shortest green path results in a new shared subpath that could be further 
improved by applying ImproveSharedSubpath recursively. 



Algorithm 3 Optimized Lensing Deformation 

1: procedure LensingDeformation(G, T) 
2: select Yj with highest weight 
3: for V e {LDU.LU,...} do 

4: = (Gi, . . . , Gy_i, decomposition V of Gy, . . . , Gy+i, . . . , G^) 

5: = (Ti, . . . , ry_i, Ty, Fy, Fy, . . . , Fy, Fy, Fy, Fy+i, . . . , F^) 

V ' 

# copies = # factors in V 

6: (G^, f ^) = SimpleDeformation(G^, F^) 

7: end for 

8: return {G^, f ^) with lowest weight 

9: end procedure 



3.3. Optimized Lensing Deformations. A single step of the optimized lensing 
deformation (Algorithm [3]) aims at improving the dominant part of the contour 
by trying various decompositions (factorizations) of its jump matrix to which, 
then, the optimized simple deformation (Algorithm [ij is applied. The contour 
parts which originate from such a lensing deformation (e.g. L, D and U in Fig. [7I 
have, however, to satisfy an additional constraint, namely condition (iv) of §[2] 
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their spatial order has to be preserved. To calculate shortest paths (line 17 in 
Algorithm [ij subject to this additional condition, we distinguish between the 
following three cases (an illustration can be found in Fig.[i3|: 






















X 




X. 



c. 



Figure 13. Illustration of the shortest path calculations done for a lensing 
deformation as in Fig. [t] Here, all three cases appear: a. All parts of the 
decomposition are still unprocessed (case 1). A shortest path pi (red) has 
to be calculated between the endpoints (magenta) of F^-. b. A path P/ for 
the contour part F| (red) to the left of F^- has already been fixed and the 
graph has been split accordingly (case 2). Now, pi must be constructed as 
the shortest path subject to the constraint that the circle pi + P/ encloses 
some point c in the split (shown in cyan). All edges are undirected, the 
arrows just indicate the directions of the paths, c. Pj (red) to the left side 
of Pi (blue) and Pr (green) to the right side of pi have already been fixed 
and the graph split accordingly (case 3). All vertices that are removed 
from the graph before calculating the shortest path pi are shown in gray. 



• Notation: T denotes the list of indices of the contour parts that are created 
by the decomposition and i denotes the index for which a shortest path is 
currently calculated. We recall that Q denotes the parts of the contour that 
have not yet been fixed in the course of Algorithm [1] 

• Case 1: T n Q = T. 

As no path belonging to T has been fixed, there are no constraints yet to 
be observed and we can simply calculate the shortest path for i. 

• Case 2: either 31 e T\Q with F/ left of F,- or 3r e T\Q with F, right of F,-. 

Without loss of generality, we assume that we have to construct a shortest 
path Ti subject to the constraint that it is to the right of an already fixed 
path P/. Now, let c be a point that is located between the left and right 
side of the split in the graph gi caused by P/. Then, the order constraint 
is identical to finding the shortest path p/ for which the circle formed by 
joining p and P/ encloses c. Lemma [1] stated at the end of this section, will 
show that Pi is actually given by 

n = {sp{gi,v\,u) ^ uw ^ sp{gi,w,v]) : uw edge in gi}, 
Pi = argmin{d(p) : p G H; ind(p + P/,c) = ±1}. 
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Hence, pi can be constructed by a minor modification of the polynomial 
algorithm O for shortest enclosing circles in embedded graphs. 

• Case 3: 31 e T\Q with T/ left of Tf and 3r e T\Q with right of T,-. 

The shortest path for F/, subject to the constraint that it is right of F/ and 
left of Tr, can only contain vertices inside the circle formed by joining P/ 
and Pr. Therefore, we construct pi as the shortest path in a smaller graph, 
in which all vertices outside of this circle have been removed. 

Lemma 1. Let q = {qi, . . . ,qn) be a path in a weighted planar graph g = {V, £), let c 
he a point considere^to he in the split along q in g[q\ and define 

n = {sp(g[^?], q\, u) + uv + s^{g[q\, v,qn) uv e E}. 

Then the shortest walk p^, in g, suhject to the constraint ind(p^ + c) = ±1, satisfies 

Pi, = argmin{d(p) : p ell; ind(p + ^,c) = ±1}. (2) 

Proof We restrict ourselves to the case ind(p + V) — because the proof for the 
other case differs just in the sign of some winding numbers. We will assume to 
the contrary that is not given by ^ and will get a contradiction. 

If there is more than one shortest walk p^, we choose the one which encloses 
the least number of vertices. If p^ i H, then there has to be a vertex v G pi, 
with p^ = I + r such that s/ = sp(g[ij], qi, v) and Sy = sp(g[ij], i^, qn) satisfy the 
following conditions 

ind(s/ + r + V/ 7^ 1/ ind(/ + Sr + '^ ,c) ^1. (3) 

We will now show that there is no such vertex v. 

Step 1. To begin with, we prove for W = p^ + ^ that 

s/ n int( W) =0, s, n int( W) = 0. (4) 

If ^ would not hold then either s/ or Sy contains a path p = (pi, . . . , pm) with 
P2/- • -/Pm-i ^ int(W) and pi,pm ^ p^- As s/ and Sr are shortest paths, such a 
subpath p has to be the shortest path from pi to pm- Consequently, the walk 

^' = PA^lrPl]+P + pAPmrqn]+X 

would satisfy 

d( WO < d( W), ind( W', c) = 1, | int( W') | < | int( W) |, 
which contradicts our choice of p^. Therefore, Q holds. 
Step 2. We now prove that 

ind( Wj, c) = 1 with Wj = / + V/ 

(5) 

ind( W[, c) = 1 with W[ = r + V. 
To this end, we consider the walk 

^We can chose any point on q for c and treat it as if it were right of q- and left of 
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which consists of the two circles 
and satisfies 

ind( W, c) = md{W\ c) = ind( Wj, c') + ind( W^, c') = 1. (6) 

If S/ n W = 0, then neither W[ nor W2 contains any vertex more than once. As g is 
a planar graph, these walks correspond to simple closed curves in the complex 
plane and therefore their winding numbers around c can only be —1, or 1. If we 
also take ^ and ([3J into account, the only possible option is 

ind( Wj,c) = 1, ind( W^,c) = 0. 

Unfortunately, s/ fi W = does not necessarily have to be satisfied. But s/ cannot 
cross / and r + ^ because of Q, which means that 

= sp(gM,^?i,i^) = sp(g[^?,T,^? + V],^?i,i;). 

Hence, we can find the following walks in g[q, 1 , + V] 

L7{ = /+ + V/ = s/ + r+ + ^, 

which are equivalent to W[ and W2 but do not contain duplicate vertices. If we 

move the paths along the splits / and + V a little bit apart, then U[ and U2 
correspond to simple connected curves, too. As we can move the split paths apart 
without crossing c or any other part of U[ or U2, these walks can only have winding 
numbers of —1, or 1 with respect to c. This means that W[ and W2 can only have 
these winding numbers even if s/ fl W 7^ 0. This proves the first relation in 
The second one follows likewise. 

Step 3. We claim that 

qn e int(W{), qn i int(W[), 

(7) 

qx e int(W[), qx i int(Wj). 

Combining ^ and ^ yields int(W) C int(Wj). Therefore all vertices in W either 
have to be in int(Wj) or in Wj. It follows that we just have to show that qn 
We consider the following closed walk 

with 

ind(WO = ind(Wj) +ind(W[) = 2. 

If qn G W{, then Si\qn,v\ = V and, consequently, contains a circle that does 
not enclose any vertex. Removing this circle from results in the walk 

W ^l\r^%\qn,qx\ 

without changing the winding number. So ind(W^') = 2, and the argument that we 
have used before to show ind(Wj, c) G { — 1, 0, 1} works for Y^" , too. Consequently, 
we get qn i y^{ - The second claim in ([7J follows likewise. 

Step 4. (see Fig. [14]) We combine the results of the previous steps to show that 
our initial assumption leads to a contradiction. As c G int(W|) fl int(W[) by (j^J, 
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the two circles W[ and W[ have a non empty intersection. Furthermore, because 
of ([tI, none of them is completely contained within the other. It follows that W[ 
and W[ have to cross each other at two or more points. One of these can be v, but 
there is actually no other vertex at which the circles could cross: s/ cannot cross 
r and, vice versa, Sr cannot cross / due to also S/ and Sr cannot cross because 
they are both shortest paths. □ 



4. Implementation Details 

4.1. The Weights. The weight de of an edge e should be an approximation of 

^||G(z)-I||d|z| 

with a suitable matrix norm. Our experiments indicate that we generally need 
fewer collocation points if we aim at minimizing all components of G — I instead of 
just focussing on its largest component, for which reason we choose the Frobenius 
(or Hilbert-Schmidt) norm. The integral is sufficiently well approximated by the 
two point trapezoidal quadrature rule (we recall that the aim of optimizing the 
weight is just preconditioning, that is, getting a particular good order of magnitude 
of the condition number). We thus take 

de = \\h-a\{\\G(h)-l\\^ + \\G{a)-l\\^) 
as the weight of an edge e with the endpoints a and h. 

4.2. The Graph. The algorithm of ^ is based on planar graphs. If the graph 
were not planar, paths could cross each other even without having any vertices 
in common and, therefore, the graph splitting described in ^ would not ensure 
that paths calculated by Algorithm [1] do not cross. We choose planar graphs built 
from rectangular grids to which a vertex in the center of each box is added that 
is connected to the vertices of the that box. Such a graph is chosen to subdivide 
a rectangle that contains all finite endpoints of F. We take this rectangle large 
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Figure 15. The rectangle to be covered by the grid is determined by the 
condition ||G — I||f > 10~^^ along F. The color coding is as in Fig.nl 

enough so that outside of it ||G — I||f is below machine precision on all arcs with 
an infinite endpoint, see Fig. [15] For numerical purposes, the jump matrix G is 
then indistinguishable from the identity matrix in the exterior of this rectangle: 
the RHP needs only to be solved in the interior. 

4.3. Contour Simplification. The algorithm described in ^returns a contour 
composed of a set of paths in the underlying graph. The collocation method of 
Olver [6], which is finally employed for the numerical solution of the RHP, would 
have to place individual Chebyshev points on each smooth (that is, linear) part of 
this piecewise linear contour. For efficiency reasons it would thus be preferable to 
have a contour with fewer breakpoints. Consequently, for each optimized path, we 
calculate a coarse piecewise linear approximation that has about the same weight. 
Quite often just a straight line connecting the endpoints of a path is already 
sufficient approximation. Fig. [i6| shows an example of this simplification process 
when applied to the final contour of Fig. [3] it cuts the number of collocation points 
by more than a factor of two while keeping the order of magnitude of the condition 
number constant. 

Conclusion 

The numerical results of this paper show that our algorithm can significantly 
reduce the condition number of RHPs. As a feature, this algorithm does not require 
any input, or knowledge, from the user other than the RHP at hand. Besides being 
thus a very convenient tool for the numerical solution of RHPs, the deformations 
automatically constructed by this algorithm might even turn out to be useful for 
determining first drafts of suitable deformations in the analytic study of RHPs. 
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